#################################################################################
#################################################################################

###                 R code to replicate simulation experiment                 ### 
###            reported in Supplementary materials (Appendix C) to            ###
###  "The Index of Emancipative Values: Measurement Model Misspecifications"  ###

#################################################################################
#################################################################################



### Required R packages ###
library(MASS)


### Simulation setup ###

# total number of simulations
nsim = 10000
# sample size
N = 100000
# number of macro-level units
n = 100
# var-cov matrix for Y1 and Y2
sigma.ind <-  matrix(c(1, 0.25, 0.25, 1), nrow = 2)
# var-cov matrix for u1 and u2
sigma.group <- matrix(c(1, 0.75, 0.75, 1), nrow = 2)
# group indicator
g <- rep(1:n, each = N/n)

# Create object to store simulation results
Results <- NULL


### Data simulation ###

for(i in 1:nsim){ 
YY <- mvrnorm(N, c(0,0), sigma.ind)
Y1 <- YY[,1]
Y2 <- YY[,2]

# country-level error
uu <- mvrnorm(100, c(0,0), sigma.group)
u1 <- rep((uu[,1]/25), each = 1000)
u2 <- rep((uu[,2]/25), each = 1000)

# Individual-level error
e1 <- rnorm(100000, 0, 1)/100
e2 <- rnorm(100000, 0, 1)/100

# country means in the absence of the macro-level method effects
mu1 <- tapply(YY[,1] + e1, g, mean)
mu2 <- tapply(YY[,2] + e2, g, mean)

#observed responses (with method effects)
y1 <- YY[,1] + e1 + u1
y2 <- YY[,2] + e2 + u2

#observed country means (with method effects)
bm1 <- tapply(y1, g, mean)
bm2 <- tapply(y2, g, mean)

Results.i <- 
cbind(cor(YY[,1], YY[,2]),        # level-1 correlation between latent preferences
      cor(y1, y2),                # level-1 correlation between observed responses
      cor(mu1, mu2),              # level-2 correlation in the absence of level-2 error
      cor(bm1, bm2), i)           # observed level-2 correlation 
Results <- rbind(Results, Results.i)
}


### Results ### 
summary(Results)

# 95%CI for the level-2 correlation in the absence of level-2 error
quantile(Results[,3], probs = c(0.025, 0.975))  # mean 0.2478

# 95%CI for the observed level-2 correlation 
quantile(Results[,4], probs = c(0.025, 0.975))  # mean 0.5547

